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ABSTRACT 



Econometricians must choose between many methods for estimating p, the 
autocorrelation coefficient, in a first order autoregressive process. This thesis examines 
the performance of four estimators in a Monte Carlo simulation. The methods 
examined arc Durbin-Watson, Beach-MacKinnon, Theil-Nagar and Prais-Winstcn. 
The autocorrelation coeficicnt, p, was varied from .2 to .9 and each method provided 
estimates of p and p, the regression coefficient, for 1000 replications. The results 
presented here are similar to those found in previous comparisons. Specifically, 
Ordinary Least Squares was found to be an efficient estimator of P when 
autocorrelation is present only to a slight degree. Of the four estimators examined, the 
performance of Theil-Nagar proved superior in estimating both p and P for small 
values of the correlation coeficicnt. Beach-MacKinnon, on the other hand, while 
containing a large bias in the estimation of p, is the more efficient estimator of P for 
large values of p. 
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I. INTRODUCTION 



A. BACKGROUND 

Autocorrelation exists in a regression model when the error terms are no longer 
independent but are correlated. In the examination of time series data autocorrelation 
is a common phenomenon and can lead to problems if Ordinary Least Squares (OLS) 
estimation procedures are used. The purpose of this thesis is to examine and compare 
four different estimates of the autocorrelation coefficient, p, the estimation of which is 
essential to the resolution of OLS deficiencies. The four estimators to be examined are 
the Durbin-Watson, Theil-Nagar, Beach-MacKinnon, and Prais-Winsten. 

B. PROBLEM STATEMENT 

In the standard regression model y=xp +e, y is a Txl vector of observations of 
a dependent variable, X is a TxK design matrix and P is a Kxl vector. The variable e 
is a Txl vector of unobservable random errors with E(e) = 0 and covariance matrix, 
E(ec') = ffMy. Thus, in the standard model, the random vector e contains elements 
which are pairwise uncorrelated with identical means and variances. In the presence of 
autocorrelation this strong assumption is violated. That is, the error terms are no 
longer independent but are correlated. The regression model becomes, 

y t = X t P + e t t= 1,2,....,T (eqn 1.1) 

where c t = p e t _j + v,, 

E(v t ) = 0, and 
E(vv') = c 2 l . 

This is known as a first order autoregressive or AR(1) process. As illustrated by 
equation 1.1, e t is expressed linearly in terms of the e t _j and another random error 
term v t . The assumption of zero mean and constant variance provides v t with all the 
nice properties of e t in the standard model. This process may occur for a variety of 
reasons, some of which arc: 

/ Omitted explanatory variables. If a correlated explanatory variable has been 
excluded from the design matrix its exclusion will be reflected in the correlation 
of the random variable e. 

2 Mispecification of the mathematical form of the model. If the wrong mathematical 
relationship is chosen the values ofe may be dependent. 
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3 Interpolations in the statistical observations. If the observational data is smoothed 
autocorrelation may result. 

4 Mispecification of the true random error. Dependence among the error terms may 
occur naturally. [Ref. l:p. 204J 

Utilizing OLS to estimate the regression coefficient, P, in the presence of an AR(1) 
process can lead to problems. Generally, there are two consequences to consider. The 
first is that the OLS estimator of the coefficients will be unbiased but will not be very 
efficient. The second consequence is that the OLS variance estimator is biased. For 
these reasons it is useful to investigate other methods to estimate P [Ref. 2:p. 439]. 

C. ESTIMATORS 

When p is known, the process is easily accounted for using Generalized Least 
Squares or Weighted Least Squares methods [Ref. 3]. However, the usual situation is 
that p is unknown and must be estimated. A number of methods have been proposed 
to estimate p and properly account for OLS deficiencies in estimating p. Chapter 2 will 
develop the four estimators mentioned above and examine the autocorrelation process. 

D. SIMULATION 

Each of the estimators considered here have the same asymptotic properties 
therefore any decision on which one to use must be based on small sample analysis and 
Monte-Carlo evidence. Therefore, a simulation will be created in which the data is 
generated according to guidelines presented in previous studies with equation 1.1 as the 
model. The actual values of p will be varied from .2 to .9. The four estimation 
techniques will then provide estimates of p and P for 1000 replications. 

E. MEASURE OF EFFECTIVENESS 

To provide an indication of which estimator performs best the mean square error 

A 4\ 

of both p and P will be estimated for each estimator. Prior results for different sets of 
estimators indicate that no one estimator will prove superior over the entire range of p 
but that one or two may out perform the others over specific intervals. 
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II. ESTIMATION 



A. GENERAL 

This chapter attempts to explore the theory behind both the first order process 
and four estimators developed to properly account for it. Three of these 
(Durbin-Watson, Thiel-Nagar, and Prais-Winsten) are categorized as estimated 
generalized least squares estimators. The fourth (Beach-MacKinnon) is a maximum 
likelihood estimator. 

B. PROCESS 

The first order process can be written as 

>' t = X t P + c t t= 1,2,...,T (eqn 2.1) 

where e, = p c t _| 4- v*, 

E(v t ) = 0, 

E(vv') = cr 2 I, 

E(v t 2 )=<r v 2 , and 
E(v t v s)=0 for s*t . 

The parameter p is generally unknown and along with |5 must be estimated. The 
statistical properties of the random error, v, listed in equation 2.1 are identical to those 
listed for e in the general linear model. The statistical properties of e under these new 
assumptions are quite different. Judge [Ref. 4:p. 438] shows that 



E(e t ) = £ p i E(v t . i ) = 0 (eqn 2.2) 

and 

E(c t 2 ) = tf e 2 = <? v 2 /( 1 -p) 2 . (eqn 2.3) 

The covariance between errors s periods apart is no longer zero and is given by 

E(e t c t _ s )= E(e t + s e t ) = (p s or 2 v )/(l-p 2 ) . (eqn 2.4) 
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The covariance matrix for e is now easily written as 



<I> = E(ec') = 



<r 2 v/0-P 2 ) 



1 P 
P 1 



T-l T-2 

p 1 1 p 1 z 



or utilizing the following convention, 



)T"1 

,T-2 



<J> = <7 v 2 T 



(eqn 2.5) 



(eqn 2.6) 



where T = 



= 1/O-P 2 ) 



1 p p 2 

p 1 p 2 

p 2 p 1 



p 

p 

p 



T-l 

T-2 

T-3 



T-l T-2 T-3 

pi i pi z pi J 



Thus, the assumptions made about the error term, e, in the standard linear model 

no longer hold for the autoregressive case. Specifically, due to autocorrelation the error 

2 • 2 

covariance matrix is no longer written as <J L l but is now <7 ^ r . 



When an attempt is made to perform a least squares fit to the data in the presence 
of an AR(1) process there arc two problems to consider. 

1 The least squares estimator p = (X'X^X'y will be unbiased but will not be very 
efficient. 

2 The least squares covariance matrix ff 2 (X'Xq‘* with "<7 2 = (y-Xb)'(y-Xb)/(T-K) 
will be a biased estimator of the variance orp. 



In the presence of positive autocorrelation Judge [Ref. 4:p. 439] notes that with 
OLS estimation the bias of the standard error of p will very likely appear as an 
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underestimate. Park and Mitchell [Ref. 5:p. 16] warn that OLS seriously 

underestimates the variance of P for p > 0.4. This understatement makes the 
estimates themselves appear much more significant than they actually are and makes 
hypothesis testing of the slope coefficients unreliable. 

C. METHODS OF ESTIMATION 

1. Generalized Least Squares Estimation 

When apriori information is available about 'F, the most convenient estimate 

A. 

for the regression coefficient, P, is obtained by applying least squares estimation 
techniques to the transformed model, 



Y* = X*p + e* (eqn 2.7) 

where Y* = PY 
X* = PX 
e* = Pe. 



The transformation matrix P is the TxT matrix 



P = 



V 1-P 2 

-p 

o 

0 



0 



0 

0 



-p 1 

0 -p 



0 0 



0 

0 

0 

1 



where P'P = T 



_ w-\ 



0 

0 

0 

0 



This method is known as the Generalized Least Squares (GLS) estimation. 

2. Estimated Generalized Least Squares 

The usual case is that p is unknown and must be estimated. Once an estimate 
✓s A . . 

for p (p) is computed one can substitute p into the P matrix and proceed with the GLS 

method outlined above. This is known as Estimated Generalized Least Squares 

(EGLS) estimation. The computational form of the alternative estimators for p 

discussed are as follows: 
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a. Durbin-Watson 
The statistic 




wlicre e t = 



(eqn 2.8) 



is often used to test for first order autoregressive errors. As the number of 
observations (T) increases it can be demonstrated that d approaches the least squares 
estimator of p or 



The Durbin-Watson statistic is provided by most least squares computer packages and 
is very easy to use. It also is an example of a two-stage estimator. That is, it first 
estimates the correlation parameter and then uses this estimate to compute the 
generalized least squares estimates for p. 
b. Theil-Nagar 

A modification of the Durbin-Watson estimator suggested by Henri 
Theil and A. L. Nagar is 



Theil and Xagar claim that this estimator is an improvement over Durbin-Watson if 
the first and second differences of the explanatory variables are small when compared 
to their corresponding ranges [Ref. 6]. Like Durbin-Watson, it also is a two-stage 
estimator. 



P = 1- (d/2) . 



(eqn 2.9) 



p- = (T 2 ( 1 -(d/2)) + K 2 ) /(T 2 - K 2 ). 



(eqn 2.10) 



c. Prais-Winsten 



A minimum sum of squares approach to estimating p yields, 




(eqn 2.1 1) 
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This estimator can be employed in both a two step and an iterative procedure. This 
paper, however, considers only the following iterative form: 

1. Set'j) = 0. 

2. Transform the variables in accordance with the transformation matrix and 
equation 2.7. 

3. Calculate the least squares estimate of P conditional on p. 

4. Calculate the estimate of p conditional on P by using equation 2.1 1. 

5. If the absolute difference in p N from the previous iteration is sufficiently small (less 
than 0.00001) stop. If not go to step 2. [Ref. 7:p. 2] 



3. Maximum Likelihood Estimation 

A maximum likelihood (ML) estimator is the value of 0 which maximizes the 
value of the likelihood function L(G). Under the assumption that Y has a multivariate 
normal distribution with mean XP and covariance matrix <r“T, the likelihood function 
is 



L(P,p,(T 2 )= C- (l/2(T 2 v )(y-Xp)'T- 1 (y-Xp) (cqn 2.12) 

where C = -(T/2)lncr 2 v + (l/2)ln (1-p 2 ) . 

The ML estimators for P, p, and cr y are those values for which, 

a L/a P = 0 , dL/dp = 0 , dL/dG 2 v = 0 . (eqn 2.13) 



Solutions to equations 2.13 are very difficult to derive. Beach and MacKinnon 
[Ref. 8:p. 54] use an ML estimator for <7 y and substitute into equation 2.12. The 
result is the concentrated likelihood function, 

L(P,p)= K-(T/2)ln((y-XP)'T' 1 (y-XP)(l-p 2 )^ T ) (eqn 2.14) 

where K = (T/2)ln(T)-(T/2) . 

They suggest maximizing L(p,p) with respect to P with p held constant and then to 
maximize with respect to p with P held constant. An algorithm to derive this ML 
estimate is 
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1. Set p = 0. 

2. Transform the variables in accordance with equation 2.7. 

An 

3. Calculate the least squares estimate of {5 conditional on p. 

4. Calculate the ML estimate of p conditional on P by solving a cubic equation of 
the untransformed residuals, (see [Ref. 8] for details) 

5. If the absolute difference in'p from the previous iteration is sufficiently small (less 
than 0.00001) stop. If not, go to step 2 [Ref. 7j. (Note: The same procedure was 
employed for iterative Prais-Winsten method except that equation 2. 1 1 was used to 
estimate p.) 

This is not a comprehensive listing of all available estimators for a first order 
process. Other estimators are listed in Judge [Ref. 4]. 
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III. COMPARISON 



A. GENERAL 

The finite sampling properties of the estimators listed here have not been derived. 
Choice of which estimator to use might be based on evidence obtained from Monte 
Carlo simulations. This chapter explains a simulation used and provides a synopsis of 
comparisons reported in the literature. 

B. PREVIOUS COMPARISONS 

There have been a number of studies of estimators for p . Each has concluded 
that OLS has serious deficiencies in the presence of autocorrelation. The majority of 
these papers have settled on two points. First, particularly in small sample sizes 
(T < 50) it is best to use estimators that consider all T observations. Rao and 
Grilitches concluded that using estimators such as Cochrane-Orcutt that ignore the 
first observation can lead to a substantial loss of efficiency [Ref. 9:p. 269j. These 
results were further substantiated by Beach and MacKinnon. In an attempt to develop 
a computationally efficient algorithm to maximize the likelihood function they 
discovered (for p= 0.6, 0.S, 0.99) significant gains in efficiency to be made by 
employing the first observation. Some of these gains are in the neighborhood of 700 
percent (Ref. 8:p. 55]. Park and Mitchell concluded that retention of this first 
observation substantially reduces the risk of collinearity as p approaches 0.9 [Ref. 5:p. 
10]. Kobayashi verified theoretically the experimental results of Park and Mitchell. By 
computing the asymptotic variances of several estimators he demonstrated that the loss 
of efficiency of the Cochrane-Orcutt method was due primarily to ignoring the first 
observation. [Ref. 10:p. 951]. 

The second point is that the Prais-Winsten solution techniques outperform many 
comparable estimators of the correlation parameter. Spitzer concluded that 
Prais-Winsten "appeared to be the best of all the two stage estimators." [Ref. ll:p. 
44]. Park and Mitchell in a later study comparing Beach-MacKinnon with the iterative 
Prais-Winsten estimator concluded that the iterative Prais-Winsten performs 
"appreciably better in estimating the autocorrelation coefficient p" [Ref. 7:p. 5[. 
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Although there were no studies found specifically comparing the four estimators 
presented here, each has demonstrated a superiority to OLS in the presence of a first 
order process. 

C. MODEL AND DATA GENERATION 

Equation 2.1 was utilized as the model with the first term in the vector e 
generated in the following fashion, 



In order to conform with previous comparisons, the data utilized in this 
experiment is identical to that used in Beach and MacKinnon [Ref. 8]. Two sample 
sizes of 20 and 50 observations were used. The untrended explanatory variable, X, was 
drawn from N(0, 0.0625) and the random error, v t , was drawn from N(0, 0.0036). 
Although autocorrelation in theory may be positive or negative, in econometric data it 
is almost’ahvays positive [Ref. l:p. 201]. For this reason p was varied from 0.2 to 0.9. 

D. VALIDATION 

The data generation program was checked to ensure the normality of e using the 
Chi Square Goodness of Fit test. The normality assumption was accepted at a 0.36S4 
level. Finally, in order to ensure each estimator performed properly the random 
portion of the model, specifically the random variable V, was removed. This allowed 
the estimators to function in a deterministic fashion. Data were then generated and 
submitted to each estimator for values of p equal 0.2, 0.6, 0.8. The results arc 
presented in Table I , illustrate that the estimators are functioning properly. 

E. SIMULATION 

For each run the values of the regression coefficients, Pq and f)j, were set to 1 
and 1. The variables X t and v f were drawn from the normal distributions discussed 
earlier. The dependent variable y t was calculated using equation 2.1. Since the 
ultimate objective was to generate residuals to send to the four estimation routines, a 
regression was then performed of y on X and residuals calculated using, 




(eqn 3.1) 



e t = y t - x t P 1 = 1.2,...,T. 



(eqn 3.2) 
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TABLE I 

ESTIMATES OF AUTOCORRELATION COEFFICIENT 



p 


DW 


TN 


PW 


BM 


.2 


.19 


.19 


.19 


.19 


.6 


.59 


.60 


.58 


.60 


.8 


.78 


oo 

O 


.77 


.80 



The values of the residuals were then sent to each estimation routine. Estimates of 
P and p were determined for values of p equal to .2, .3, .4, .5, .6, .7, .8, and .9. Each 
estimate was replicated 1000 times for the sample sizes of 20 and 50. 

F. MEASURES OF EFFECTIVENESS 

In order to compare the performances of the estimators, two MOE's were used. 
The mean square error (MSF) of p was estimated for each estimator. This represents 
the expected squared error made in estimating p. The following computational form of 
MSE was used, 



/OCO 

E (P - Pj)/' 1000 i= 1 ,..., 1000. 

/-/ 



(cqn 3.3) 



The successive values of MSE of p were then plotted against the actual p to 
examine performance over the range of p. 

The second MOE examined the relative^efficicncies of the regression coefficient a^ 
defined in [Ref. 7:p. 7]. A ratio of MSE of P for a particular estimate to the MSE of P 
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for the OLS estimate allows the examination of the relative gains in using particular 
techniques over OLS. Since the proper estimation of P is paramount the efficiency of P 
is predetermined to be the most important MOE. 
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IV. RESULTS AND CONCLUSIONS 



A. GENERAL 

The major emphasis of this thesis was to examine the performance of four 
estimators of the autocorrelation coefficient, p, for a first order autoregressive process. 
The estimators examined were Durbin-Watson, Theil-Nagar, Prais-Winsten, and 
Beach-MacKinnon. 

A Monte-Carlo simulation was performed for the following values of p: .2, .3, .4, 
.5, .6, .7, .8, and .9. Each run was replicated 1000 times for sample sizes of 20 and 50. 
The results are recorded in Table 11 . Irrespective of sample size, each of the methods 
underestimate the true value of p but as the number of observations is increased from 
20 to 50 the bias reduces. As was expected, no one estimator uniformly outperforms 
the others. In both sample sizes, the two stage estimators (Durbin-Watson and 
Theil-Nagar) achieve better results for small p. As the value of p increases, the 
iterative methods (Prais-Winsten and Beach-MacKinnon) perform best. With T=20 
this transition occurs at p= .6 while at 50 observations it occurs earlier at p = .4. 

The discussion of the results will be divided into two sections. The measures of 
effectiveness, as defined in Chapter 3 will first be applied to the simulation results for 
T = 20. This will be followed by an identical approach when the sample size is 
increased to 50. 

B. SAMPLE SIZE 20 

Since performance of an estimator is roughly indicated by its mean and variance, 
mean square error (MSE) of each p* over the entire sample size was estimated. The 
results of these calculations are presented in Table III along with a plot of MSE of p^ 
versus actual values of p in Figure 4.1 . They again indicate that the Theil-Nagar and 
Durbin-Watson estimators are better for smaller values of p (p<.6) and as p increases 
the Prais-Winsten p emerges as the best. On the basis of Figure 4.1 alone, 
Bcach-MacKinnon's performance is clearly inferior. However, in examining the 
efficiency of each estimator in Table IV, Beach-MacKinnon proves to be the most 
efficient in estimating |3 over the widest range of p. The tie in Figure 4.1 between 
Theil-Nagar and Durbin-Watson is resolved in Table IV with Theil-Nagar proving to 
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TABLE 11 

ESTIMATES OF AUTOCORRELATION COEFFICIENT 



Sample Size 20 

P 

.2 

.3 

.4 

.5 

.6 

.7 

.8 

.9 

Sample Size 50 

P 

.2 

.3 

.4 

.5 

.6 

.7 

.8 

.9 



DW 


TN 


.158 


.162 


.234 


.239 


.310 


.316 


.385 


.392 


.460 


.467 


.533 


.542 


.603 


.613 


.667 


.681 


DW 


TN 


.173 


.161 


.279 


.253 


.360 


.340 


.451 


.432 


.544 


.523 


.630 


.610 


.726 


.700 


.812 


.796 



PW 


BM 


.113 


.107 


.205 


.193 


.296 


.279 


.387 


.365 


.478 


.450 


.567 


.535 


.655 


.617 


.741 


.697 


PW 


BM 


.170 


.178 


.270 


.270 


.360 


.359 


.463 


.453 


.559 


.547 


.651 


.640 


.747 


.731 


.839 


.820 
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be uniformly more efficient than Durbin-Watson. Table IV also demonstrates that for 
p = .2 OLS is at least as efficient as three of the four estimators. 

C. SAMPLE SIZE 50 

The results of the MSE calculations for T=50 arc recorded in Table V along 
with a plot of MSE of p versus the actual values of p in Figure 4.2. The 
Durbin-Watson and Theil-Nagar estimators again perform the best for smaller values 
of p (p < .4) and as p increases the Bcach-McKinnon and Prais-Winstcn estimators of 
p contain the smallest MSE. 

Once again even though the Prais-Winsten p has a smaller MSE than 
Bcach-MacKinnon, Table VI illustrates that Bcach-McKinnon is a uniformly more 
efficient estimator of the slope coefficient. For the smaller values of p (p < .4) 
Theil-Nagar is more efficient than Durbin-Watson. Table VI also illustrates that OLS 
is at least as efficient as any of the other estimators when p is small. 

D. SUMMARY 

As was found in previous studies when autocorrelation is present only to a slight 
degree (p<.2) the OLS estimator provides an efficient estimate for the regression 
coefficient, p. As the process becomes more significant however, all the estimators 
outperform the OLS solution. In both sample sizes the performance of Theil-Nagar 
and Durbin-Watson are nearly identical with respect to the MSE of p. However, when 
efficiency of the slope coefficient estimate is examined, Theil-Nagar proves to be the 
better 2 stage estimator. Park and Mitchell [Ref. 7:p. 4] found that Prais-Winstcn 
performs better in estimating p. The results presented here tend to dispute that 
finding. For while Prais-Winstcn has a uniformly smaller MSE of p, 
Bcach-MacKinnon provides the most efficient estimator of P . Spitzcr, on the other 
hand (Ref. ll:p. 44], which ranked two stage estimators as being the best for values of 
p between .2 and .5, mirrors the results produced here. Apriori knowledge of the 
neighborhood of p will be helpful in selecting the appropriate estimation method. For 
both sample sizes Theil-Nagar appears to be the best for small values of p. 
Bcach-MacKinnon, while containing a larger bias for p than docs Prais-Winsten, is a 
much more efficient estimator of the slope coefficient for larger values of p. 
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Figure 4.1 Estimated mean square error of p vs. p (sample size = 20). 



TABLE 111 

DATA PRESENTED IN FIGURE 4.1 
Sample Si/e 20 



p 


MSEDW 


MSETN 


MSEPW 


MSEBM 


.2 


.7494 


.7485 


.8557 


.8594 


.3 


.6268 


.6244 


.7004 


.7115 


.4 


.5156 


.5124 


.562 1 


.5787 


.5 


.4159 


.4126 


.4400 


.4604 


.6 


.3278 


.3250 


.3342 


.3566 


.7 


.2519 


.2495 


.2433 


.2662 


.8 


.1891 


.1860 


.1698 


.1916 


.9 


.1407 


.1357 


.1141 


.1331 
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TABLE IV 

EFFICIENCY OF REGRESSION COEFFICIENT ESTIMATES 
Sample Si/e 20 



l> 


MSEp (1)W) 


MSEp (TN) 


MSEp (I»\V) 


MSEp (BM) 




MSEp (OLS) 


MSEp (OLS) 


MSEp (OLS) 


MSEP (OLS) 


7 

« i- 


1 .004 


.9794 


1.035 


1. 04 1 


.3 


.9228 


.8967 


.9442 


.9515 


.4 


.8218 


.7929 


.8325 


.8342 


.5 


.7082 


.6751 


.7024 


.6959 


.6 


.5864 


.54S4 


.5652 


.5515 


.7 


.4610 


.4207 


.4329 


.4135 


.8 


.3359 


.3020 


.3093 


.2S70 


.<) 


.2253 


.2087 


.2077 


.1892 
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Figure 4.2 Estimated mean square error of p vs. p (sample size = 50). 



TABLE V 

DATA PRESENTED IN FIGURE 4.2 
Sample Size 50 



p 


MSEDW 


MSETN 


MSEPW 


MSEBM 


.2 


.7010 


.6766 


.7055 


.7065 


.3 


.5500 


.5399 


.5653 


.5578 


.4 


.4383 


.4196 


.4396 


.5578 


.5 


.3298 


.3151 


.3056 


.3156 


.6 


.2358 


.2262 


.2116 


.2215 


.7 


.1578 


.1526 


.1357 


.1449 


.8 


.0906 


.0942 


.0773 


.0851 


.9 


.0500 


.0509 


.0360 


.0417 
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TAB LIZ VI 

EFFICIENCY OF REGRESSION COEFFICIENT ESTIMATES 
Sample Size 50 



p 


MSFp (DW) 


MSFP (TN) 


MSFP (PW) 


MSEp (DM) 




MSFp (OLS) 


MSFP (OFS) 


MSFp (OFS) 


MSFp (OFS) 


.2 


1.073 


1.041 


1.046 


1.058 


.3 


.9985 


.9482 


.9714 


.9562 


.4 


.8850 


.8255 


.8635 


.8255 


.5 


.7452 


.6859 


.7282 


.6825 


.6 


.5920 


.5420 


.5870 


.5406 


.7 


.4366 


.4020 


.4453 


.4020 


.8 


.2889 


.2690 


.3067 


.2700 


.9 


.1589 


.1505 


. 1 73S 


.1505 
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APPENDIX 
PROGRAM LISTINGS 



This appendix contains listings of the programs utilized in the analysis performed 
herein. All of the functions are written in FORTRAN and contain the necessary 
documentation. The Monte Carlo simulation was performed using the Advanced 
Simulation and Statistics Package [Ref. 12] developed by P. A. Lewis. Since the 
package only allows for the simultaneous comparision of 3 estimators, 2 functions were 
developed for each sample size. The first, SIMS generates estimates for 
Durbin- Watson, Thcil-Nagar, and Prais-Winsten for a sample size of 20. SI MSA 
meanwhile, generates estimates for Beach-MacKinnon for the identical sample size. 
Routines for Durbin-Watson and Theil-Nagar were included in SI MSA to ensure the 
results were comparable to SIMS. SI MSB and SIMSC perform in a similar fashion for 
sample size of 50 and therefore were not included. The Advanced Simulation and 
Statistics Package computes the mean square error of "p" for each estimator 
automatically. The mean square error for the p estimates was accomplished by the 
MSEB function. 



SIMS 

DIMENSION EHAT(20) 

COMMON /MYDATA/ K,T,ANS,Y1,X 

COMMON /DATA1/ IX1A.RH0 

REALM Y( 5000 ) , YMIN , YMAX , PMEAN(3) 

CHARACTERS T1,T2,T3 

INTEGER N, M,NE(8),L,D,RG,SEI,SVS, NEST, NSR.IXl, 1X2, 1X3 

EXTERNAL DATGEN, DURWAT, BEAMAC, PRAWIN, LSEB, DCALC , TRANSF 

EXTERNAL LN0RM , S IMTBD , GMPRD 

NR=20 

T=20 

K=2 



C 

C 
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c 

OPEN (UNIT=19 , FI LE=‘ MONICA 1 ) 

C OPEN ( UN IT=21 , FI LE= ' MARGE 1 ) 

C OPEN (UNIT=51 , FI LE =I AMBROSE 1 ) 

C OPEN (UNIT=41 , FI LE= ' DAT 2 1 ) 

OPEN (UN IT=6 1 , FI LE= 1 DAT3 1 ) 

READ (19,*) ANS 

10 READ( 19 ,* , END=999) N,M,L,D,RG,SEI ,SVS,NEST,NSR 
READ(19,*)YMIN,YMAX 
READ( 19,*) (NE(I) ,1=1 ,L) 

READ( 19 , 120) 1X1 , 1X2 , 1X3 
120 FORMAT (I5,1X,I5,1X,I5) 

READ (19,115) T1 
115 F0RMAT(A80) 

READ( 19 ,115) T2 
READ ( 19 , 11 5)T3 
READ( 19 ,*) (PMEAN(I) ,1=1,3) 

READ( 19 ,*) RHO 
READ (19,61)IX 1 A 
61 FORMAT(I5) 

C 

C CALL FOR SIMTBD 

C 

CALL SIMTBD ( 1X1 , 1X2 , 1X3 , Y ,N ,M,NE , L,D ,NSR, RG ,SEI ,SVS, 

*YMIN , YMAX , NEST , DATGEN , DURWAT , T1 , DATGEN , BEAMAC ,T2 , DATGEN , PRAWIN ,T3 , 
*PMEAN) 

GO TO 10 

999 WRITE(6,*)'END OF DATA INPUT 1 

STOP 
END 

q **************************************************************** 

q ************* q/\TA generation subroutine ********************** 

0 ***************************************************************** 

SUBROUTINE DATGEN (1X1 ,EHAT,NR) 

DIMENSION BHAT(2) ,YSTAR(20) ,R2(20) ,U(20) , 



2S 



*E(20) ,YHAT(20) ,EHAT(20) ,XSTAR(20,2) 

*,Y1(20) ,X(20,2) , V(20) 

COMMON /MYDATA/ K,T,ANS,Y1,X 
COMMON /DATA1/ IX1A.RHO 
INTEGER 1X1 , IX1A,NR 
C 
C 

C GENERATE THE RANDOM ERROR 

C 

CALL SNOR ( 1X1 , U , NR ,1,0) 

C 

C ADJUST THE VARIANCE OF R. E. IAW BEACH AND MACKINNON( 1978) 

DO 38 1=1 ,T 

V(I)=U(I)\06 
38 CONTINUE 

C 

C GENERATE THE ERROR FOR THE STAND LINEAR MODEL 

C 

E( 1 )=V( 1 )/( 1-(RH0**2) )**0. 5 
DO 31 J=2,T 

E(J)=RHO*E( J-1)+V(J) 

31 CONTINUE 
C 

C 

C GENERATE THE EXPLANATORY VARIABLES IAW RAO AND GRILITCHES (1969) 
C 

DO 32 1=1,20 

X(I,1)=1 

32 CONTINUE 

C CHANGE 1X1 IN ORDER TO AVIOD COLLINEARITY 

C I X 1 A= I X 1 + 1 9 

CALL SN0R(IX1A,R2,NR, 1 ,0) 

DO 33 J=1 ,20 

X(J ,2)=R2( J)*. 25 

33 CONTINUE 
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THE TRUE BETA EQUALS 1,1 



C 
C 
C 
C 
C 

C GENERATE THE INDEPENDENT VARIABLE 
C 

DO 35 1=1,20 

Y1(I)=(X(I,1)+X(I,2))+E(I) 

35 CONTINUE 

C 

C GENERATE THE LEAST SQUARES ESTIMATOR FOR BETA 

C 

CALL LSEB(X,Y1,BHAT) 

C PRINT LSEB TO A FILE 

IF(ANS .EQ. 2) WRITE(61,201)BHAT 
201 FORMAT( Fll. 8,2X , Fll. 8) 

C 

C GENERATE YHAT 

C 

DO 100 1=1,20 

YHAT ( I )=X( I ,1)*BHAT(1)+X(I ,2)*BHAT(2) 

100 CONTINUE 

C 

c 

C GENERATE EHAT 

C 

DO 50 1=1,20 

EHAT ( I )=Y 1(1) -YHAT ( I ) 

50 CONTINUE 

C 

C 

RETURN 

END 

q ********************* DURBIN WATSON *************************** 
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C THIS FUNCTION COMPUTES THE DURBIN-WATSON ESTIMATE OF RHO 
REAL FUNCTION DURWAT (EHAT.NR'.WI) 

DIMENSION EHAT(20) ,X(20 ,2) , Yl( 20 ) , XSTAR1 ( 20,2) , YSTAR1(20) ,BHAT1(2) 
COMMON /MYDATA/ K,T,ANS,Y1,X 
CALL DCALC (EHAT,T,D) 

DURWAT=l-D/2 

C 

CALL TRANSF(X,Y1 , DURWAT ,XSTAR1,YSTAR1) 

CALL LSEB (XSTAR1 ,YSTAR1 ,BHAT1) 

IF (ANS . EQ. 1 ) WRITE(21, 701) BHAT1 
701 F0RMAT(F11. 8,2X,F11. 8) 

C 

C 

C 

RETURN 

END 

C 

c 

£ ** * * * ** * * * * ** ** ** ** * THEIL NAGAR * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

C THIS FUNCTION COMPUTES THE THEIL-NAGAR ESTIMATE OF RHO 

C REAL FUNCTION THENAG (EHAT.NR.WI) 

DIMENSION EHAT(20) , YSTAR2(20) ,XSTAR2(20 ,2) ,BHAT2(2) 

*,Y1(20) , X( 20 , 2 ) 

COMMON /MYDATA/ K,T,ANS,Y1,X 
CALL DCALC (EHAT,T,D) 

THENAG=((T**2)*( l-D/2)+K**2)/(T**2-K**2) 

CALL TRANSF(X ,Y1 , THENAG, XSTAR2 ,YSTAR2) 

CALL LSEB ( XSTAR2 , YSTAR2 , BHAT 2 ) 

IF (ANS .EQ. 1 ) WRITE(31 ,801) BHAT 2 
801 F0RMAT(F11. 8,2X,F11. 8) 

RETURN 

END 

q *********************** p VJINSTEN ************************* 

C THIS FUNCTION COMPUTES THE PRAIS-WINSTEN ESTIMATE OF RHO 
REAL FUNCTION PRAWIN(EHAT,NR,WI) 
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DIMENSION EHAT3(20) ,YHAT3(20) ,YSTAR3(20) ,BHAT3(2) , 
*EHAT(20) ,XSTAR3(20 ,2) 

* ,Y1(20) ,X(20 ,2) 

COMMON /MYDATA/ K,T,ANS,Y1,X 
N=0 

RH03=0 
98 N=N+I 

CALL TRANSF (X , Y1 , RH03 .XSTAR3 , YSTAR3) 

CALL LSEB (XSTAR3 ,YSTAR3,BHAT3) 

C GENERATE YHAT3 

DO 83 1=1,20 

YHAT3( I )=X( I , 1)*BHAT3( 1 )+X( I ,2)*BHAT3(2) 

83 CONTINUE 

DO 4 1=1, T 

EHAT3(I)=Y1(I)-YHAT3(I) 

4 CONTINUE 
C 

RH0NUM=0 
RHODEN=0 
DO 5 1=2, T 

RHONUM=RHONUM+( EHAT3( I )*EHAT3( 1-1 ) ) 

5 CONTINUE 
C 

DO 6 1=2, T-l 

RHODEN=RHODEN+( EHAT3( I )**2) 

6 CONTINUE 

PRAWIN=RHONUM/RHODEN 

C CHECK FOR PRAWIN WHICH ARE OUT OF BOUNDS 

I F ( P RAW I N . GE. 1)THEN 
PRAWIN=0. 99999 
ELSE IF (PRAWIN. LE.-1)THEN 
PRAWIN=-0. 99999 
END IF 

C COMPARISON OF RH03 AND PRAWIN IF DIFF . LT. 0.0001 THEN END 
IF(ABS( RH03- PRAWIN). GT. . 0001)THEN 



bl 



RH03=PRAWIN 
GO TO 98 
ELSE 

PRAWIN=PRAWIN 
END IF 

C IF (ANS .EQ. 1 ) WRITE(41, 901) BHAT3 
C01 FORMAT( Fll. 8,2X , Fll. 8) 

RETURN 

END 

C 

C THE FOLLOWING SUBROUTINES AID IN THE COMPUTATION OF THE FOUR 
C ESTIMATORS OF RHO. 

q ************* SUBROUTINE LSEB ******************************* 

C SUBROUTINE LSEB WILL COMPUTE THE LSE OF B 
C 

SUBROUTINE LSEB(X ,Y1 ,BHAT) 

DIMENSION BHAT(2) ,Y1(20),X(20,2) ,XTRNSP(2,20) ,XI(2,2),H(2,20), 
*XPRIX(2,2) 

C X TRANSPOSE 

DO 40 1=1,20 
__ DO 41 J=1 ,2 

XTRNSP(J , I)=X( I ,J) 

41 CONTINUE 

40 CONTINUE 

C MULTIPLY X TRANSPOSE AND X 

CALL GMPRD(XTRNSP,X,XPRIX,2,20,2) 

C CALCULATE INVERSE OF X PRIME X 

DETR=1/(XPRIX( 1,1)*XPRIX(2,2)-XPRIX(1,2)*XPRIX(2,1)) 

XI( 1 ,1)=DETR*XPRIX(2,2) 

XI(1,2)=DETR*(-XPRIX(1,2)) 

XI(2,1)=DETR*(-XPRIX(2,1)) 

XI(2,2)=DETR*XPRIX( 1,1) 

C MULTIPLY INVERSE AND TRANSPOSE 

CALL GMPRD( X I ,XTRNSP,H,2,2,20) 

DO 99 1=1,2 



33 



BHAT( I)=H( I ,1)*Y1(1)+H(I,2)*Y1(2)+H(I,3)*Y1(3) 
*+H(I,4)*Yl(4)+H(I,5)*Yl(5) 

*+H( I,6)*Y1(6)+H(I,7)*Y1(7)+H( I ,8)*Y1(8)+H( I,9)*Y1(9) 

*+H( I , 10)*Y1( 10)+ 

*H(I,11)*Y1(11)+H(I,12)*Y1(12)+H(I,13)*Y1(13)+H(I,14)*Y1(14) 
*+H( I , 15)*Y1( 15)+ 

*H(I,16)*Y1(16)+H(I,17)*Y1(17)+H(I,18)*Y1(18)+H(I,19)*Y1(19)+ 

*H(I,20)*Y1(20) 

99 CONTINUE 

RETURN 
END 

q ************* SUBROUTINE DCALC ***************************** 

C SUBROUTINE DCALC WILL COMPUTE THE DURBIN STATISTIC D 

C 

SUBROUTINE DCALC(EHAT,T,D) 

DIMENSION D1(20),D2(20), EHAT(20) 

DNUM=0 
DDEN=0 
DO 1 1=2, T 

D1(I-1)=(EHAT(I)-EHAT(I-1))**2 
DNUM=DNUM+D1( 1-1) 

1 CONTINUE 

DO 2 J=1,T 

D2(J)=EHAT(J)**2 

DDEN=DDEN+D2(J) 

2 CONTINUE 

D=DNUM/DDEN 

RETURN 

END 

C 

q ***************** SUBROUTINE TRANSF ********************** 

c 

C SUBROUTINE TRANSF IS DESIGNED TO TRANSFORM THE X'S AND Y'S 
C ACCORDING TO THE LEAST SQUARES RULE. 

SUBROUTINE TRANSF( X , Y1 , RHOHAT , XSTAR , YSTAR) 
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DIMENSION Yl(20) ,YSTAR(20) ,X(20,2) , XSTAR(20,2) 

K=2 

T=20 

C Y TRANSFORM 

YSTAR( 1)=(( 1-(RHOHAT**2))**0. 5)*Y1(1) 

DO 7 1=2,20 

YSTAR(I)=Y1(I)-(RH0HAT*Y1(I-1)) 

7 CONTINUE 

C X TRANSFORM 

XSTAR( 1 , 1)=( 1-(RHOHAT**2))**0. 5 . 

DO 9 J=2,K 

XSTAR(1,J)=((1-(RHOHAT**2))**0. 5)*X(1,J) 

9 CONTINUE 

DO 11 L=2,T 

XSTAR(L,l)=l-RHOHAT 

11 CONTINUE 

DO 12 1=2, T 
DO 13 J=2, K 

XSTAR( I , J)=X( I , J)-RHOHAT*X( 1-1 , J) 

13 CONTINUE 

12 CONTINUE 

RETURN 

END 
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S I MSA 

C THE PURPOSE OF THIS PROGRAM IS TO RUN COMPUTE THE FOLLOWING 

C ESTIMATORS (DW TN BM) FOR A SAMPLE SIZE OF 20 

DIMENSION EHAT(20) 

COMMON /MYDATA/ K,T,ANS,Y1,X 

COMMON /DATA1/ IX1A.RH0 

REAL*4 Y ( 5000 ) ,YMIN,YMAX,PMEAN(3) 

CHARACTER*80 T1,T2,T3 

INTEGER N, M,NE(8),L,D,RG,SEI,SVS, NEST, NSR,IX1, 1X2, 1X3 
EXTERNAL DATGEN, DURWAT, THENAG, BEAMAC, LSEB, DCALC , TRANS F 
EXTERNAL LNORM,SIMTBD,GMPRD 
NR=20 
T=20 
K=2 
C 
C 

c 

OPEN ( UN IT=19, FI LE=' MONICA 1 ) 

OPEN ( UN IT=51, FI LE =I AMBROSE 1 ) 

READ (19,*) ANS 

10 READ( 19 ,* , END=999) N,M,L,D,RG,SEI,SVS,NEST,NSR 
READ(19,*)YMIN,YMAX 
READ( 19 ,*) ( NE( I ) , 1=1 , L) 

WRITE (22,105) (NE( I) , 1=1 , L) 

105 F0RMAT(8I4) 

READ(19,120) 1X1 ,1X2 ,1X3 
120 FORMAT( 15, IX, 15, IX, 15) 

READ (19,115) T1 
115 F0RMAT(A80) 

READ(19,115) T2 
READ ( 19 ,115 )T3 
READ(19,*) (PMEAN( I) , 1=1,3) 

READ( 19 ,*) RHO 
READ( 19 , 6 1 ) I X 1A 
61 FORMAT( 15) 
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CALL FOR SIMTBD 



C 
C 
C 

CALL SIMTBD (IX1,IX2,IX3,Y,N,M,NE,L,D,NSR,RG,SEI,SVS, 

*YMIN, YMAX , N EST , DATGEN , DURWAT , T1 , DATGEN , THENAG , T2 , DATGEN , BEAMAC , T3 , 
*PMEAN) 

GO TO 10 

999 WRITE(6 ,*) 1 END OF DATA INPUT' 

STOP 

END 

Q **************************************************************** 

q ************* QA7A GENERATION SUBROUTINE ********************** 

0 ***************************************************************** 

SUBROUTINE DATGEN (IX1,EHAT,NR) 

DIMENSION BHAT(2) ,YSTAR(20) ,R2(20),U(20), 

*E(20) ,YHAT(20) ,EHAT(20) ,XSTAR(20 ,2) 

* ,Y1(20) ,X(20 ,2) 

COMMON /MYDATA/ K,T,ANS,Y1,X 
COMMON /DATA1/ IX1A.RHO 
INTEGER 1X1 , IX1A.NR 
C _ 

C 

C GENERATE THE RANDOM ERROR 

C 

CALL SNOR ( 1X1 , U ,NR , 1 ,0) 

C 

C 

C GENERATE THE ERROR FOR THE STAND LINEAR MODEL 

C 

E(l)=U(l)/( l-( RHO**2) )**0. 5 
DO 31 J=2,20 

E(J)=RHO*E(J-l)+U( J) 

31 CONTINUE 

C 

C 
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C GENERATE THE EXPLANATORY VARIABLES I AW RAO AND GRILITCHES (1969) 
C 

DO 32 1=1,20 

X(I,1)=1 

32 CONTINUE 

C CHANGE 1X1 IN ORDER TO AVIOD COLLINEARITY 

C IX1A=IX1+19 

CALL SNOR( IX1A,R2 ,NR, 1 ,0) 

DO 33 J=1 ,20 

X(J ,2)=R2(J)*. 25 

33 CONTINUE 
C 

C 

C THE TRUE BETA EQUALS 1,1 

C 

C 

C GENERATE THE INDEPENDENT VARIABLE 

C 

DO 35 1=1,20 

Yl( I )=( X( I , 1)+X( I ,2))+E( I) 

35 CONTINUE - 

C 

C 

C GENERATE YHAT 

C 

CALL LSEB(X,Y1,BHAT) 

C BHAT ( 1 ) = 1 . 3 

C BHAT(2)=1. 1 

DO 100 1=1,20 

YHAT ( I )=X ( I , 1)*BHAT(1)+X(I ,2)*BHAT(2) 

100 CONTINUE 

C 

C 

C GENERATE EHAT 

C 



3S 



50 

C 

C 



C 

C 



C 

c 

c 



c 

c 

c 

c 



c 

c 



DO 50 1=1,20 

EHAT(I)=Y1(I)-YHAT(I) 

CONTINUE 



RETURN 

END 

* * * * * * * * * ** * * rt * * rt rtrt * * DURBIN WATSON *************************** 
REAL FUNCTION DURWAT (EHAT,NR,WI) 

DIMENSION EHAT(20) ,X(20,2),Y1(20) ,XSTAR1(20 ,2) , YSTAR1(20) ,BHAT1(2) 
COMMON /MYDATA/ K,T,ANS,Y1,X 
CALL DCALC (EHAT,T,D) 

DURWAT=l-D/2 

CALL TRANSF(X,Y1, DURWAT, XSTAR1.YSTAR1) 

CALL LSEB (XSTAR1 , YSTAR1 ,BHAT1) 



RETURN 

END 



****** ************** JHEIL NAGAR ****************************** 

REAL FUNCTION THENAG (EHAT,NR,WI) 

DIMENSION EHAT(20) ,YSTAR2(20) ,XSTAR2(20,2) ,BHAT2(2) 
*,Y1(20),X(20,2) 

COMMON /MYDATA/ K,T,ANS,Y1,X 
CALL DCALC (EHAT,T,D) 

THENAG=( (T**2)*(l-D/2)+K**2)/(T**2~K**2) 

RETURN 

END 



BEACH MACKINNON 



*********************** 
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************************* 



c 

REAL FUNCTION BEAMAC( EHAT , NR , WI ) 

DIMENSION EHAT4( 20 ) , YHAT4( 20 ) ,YSTAR4(20) ,BHAT4(2) , 
*Y1(20) ,EHAT(20) ,X(20,2) , XSTAR4(20 , 2) 

COMMON /MYDATA/ K,T,ANS,Y1,X 
N=0 

RH04=0 
98 N=N+1 

CALL TRANSF (X , Y1 , RH04 ,XSTAR4 , YSTAR4) 

CALL LSEB ( XSTAR4 , YSTAR4 , BHAT4) . 

C BHAT4( 1)=1. 0 

C BHAT4(2)=1. 0 

C GENERATE YHAT4 

DO 83 1=1,20 

YHAT4( I )=X( I , 1)*BHAT4( 1)+X( I , 2)*BHAT4(2) 

83 CONTINUE 

DO 4 1=1, T 

EHAT4( I )=Y1( I)-YHAT4( I ) 

4 CONTINUE 

SUM3=0 
SUM2=0 
SUM1=0 
DO 71 1=2, T 

SUMl=SUMl+( EHAT4( I ) *EHAT4( 1-1 ) ) 

71 CONTINUE 
C 

DO 72 1=2, T 

SUM2=SUM2+(EHAT4(I-1)**2) 

72 CONTINUE 
C 

DO 73 1=2, T 

SUM3=SUM3+( EHAT4( I )**2) 

73 CONTINUE 
C 

DEN0M=(T-1)*(SUM2-(EHAT4(1)**2)) 
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c 

A=(-(T-2)*SUM1)/DEN0M 

C 

B=(((T-1)*(EHAT4( 1)**2))-(T*SUM2)-SUM3)/DEN0M 
C 

C=(T*SUM1)/DEN0M 

C 

SMALP=B-( ( A**2)/3) 

C 

SMALQ=C-((A*B)/3)+( (2*(A**3))/27) 

C 

THETA=ACOS((SMALQ*(27**. 5))/(2*SMALP*((-SMALP)**0. 5))) 

C 

C BEAMAC IS THE ITERATIVE RHO FOR THIS PROCEEDURE 

BEAMAC=(-2*( (-SMALP/3)**0. 5))*COS((THETA/3)+(3. 1412/3 )) -(A/3 ) 
C CHECK FOR BEAMAC WHICH ARE OUT OF BOUNDS 
I F( BEAMAC. GE. 1)THEN 
BEAMAC=0. 99999 
ELSE IF (BEAMAC. LE. -1)THEN 
BEAMAC=-0. 99999 

END IF — 

C COMPARISON OF RH04 AND BEAMAC IF DIFF . LT. O.OOOl THEN END 
IF(ABS(RH04-BEAMAC). GT. . 0001)THEN 
RH04=BEAMAC 
GO TO 98 
ELSE 

BEAMAC=BEAMAC 
. END IF 

IF (ANS .EQ. 2) WRITE (51,901) BEAMAC 
901 FORMAT( F15.ll) 

RETURN 

END 

C 

C THE FOLLOWING SUBROUTINES AID IN THE COMPUTATION OF THE FOUR 
C ESTIMATORS OF RHO. 
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q * ** ***** ***** SUBROUTINE LSEB ******************************* 
C SUBROUTINE LSEB WILL COMPUTE THE LSE OF B 
C 

SUBROUTINE LSEB(X ,Y1 ,BHAT) 

DIMENSION BHAT(2) ,Y1(20),X(20,2) , XTRN SP ( 2 ,20),XI(2,2),H(2,20), 
*XPRIX(2 ,2) 

C X TRANSPOSE 

DO 40 1=1,20 
DO 41 J=1 ,2 

XTRNSP(J ,I)=X(I 
41 CONTINUE 

40 CONTINUE 

C MULTIPLY X TRANSPOSE AND X 

CALL GMPRD(XTRNSP,X ,XPRIX ,2 ,20 ,2) 

C CALCULATE INVERSE OF X PRIME X 

DETR=1/(XPRIX( 1 , 1)*XPRIX(2 ,2)-XPRIX( 1 ,2)*XPRIX(2 , 1) ) 

X I ( 1 , 1 )=DETR*XPRIX(2 , 2) 

XI( 1 ,2)=DETR*(-XPRIX( 1,2)) 

XI(2,1)=DETR*(-XPRIX(2,1)) 

XI(2 ,2)=DETR*XPRIX( 1,1) 

C MULTIPLY INVERSE AND TRANSPOSE - 

CALL GMPRD (XI ,XTRNSP ,H,2 ,2,20) 

DO 99 1=1,2 

BHAT ( I )=H ( I , 1)*Y1( 1)+H( I , 2)*Y1(2)+H( 1 ,3)*Y1(3) 
*+H(I,4)*Yl(4)+H(I,5)*Yl(5) 

*+H( I ,6)*Y1(6)+H(I,7)*Y1(7)+H(I,8)*Y1(8)+H(I,9)*Y1(9) 
*+H(I,10)*Yl(10)+ 

*H(I,11)*Y1(11)+H(I,12)*Y1(12)+H(I,13)*Y1(13)+H(I,14)*Y1( 14) 
*+H(I,15)*Yl(15)+ 

*H(I,16)*Y1(16)+H(I,17)*Y1(17)+H(I,18)*Y1(18)+H(I,19)*Y1(19)+ 
*H( I ,20)*Y1(20) 

99 CONTINUE 

RETURN 
END 

q ************* SUBROUTINE DCALC ***************************** 
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C SUBROUTINE DCALC WILL COMPUTE THE DURBIN STATISTIC D 
C 

SUBROUTINE DCALC(EHAT,T,D) 

DIMENSION Dl(20) ,D2(20) , EHAT(20) 

DNUM=0 
DDEN=0 
DO 1 1=2, T 

D1(I-1)=(EHAT(I)-EHAT(I-1))**2 

DNUM=DNUM+D1(I-1) 

1 CONTINUE 

DO 2 J=1,T 

D2(J)=EHAT(J)**2 

DDEN=DDEN+D2(J) 

2 CONTINUE 

D=DNUM/DDEN 

RETURN 

END 

C 

q ***************** SUBROUTINE TRANS F ********************** 

c 

C SUBROUTINE TRANS F IS DESIGNED TO TRANSFORM THE X'S AND Y'S 
C ACCORDING TO THE LEAST SQUARES RULE. 

SUBROUTINE TRANSF( X , Y1 , RHOHAT , XSTAR , YSTAR) 

DIMENSION Y 1 ( 20 ) ,YSTAR(20) ,X(20,2) ,XSTAR(20,2) 

K=2 

T=20 

C Y TRANSFORM 

YSTAR( 1)=(( 1-(RH0HAT**2) )**0. 5)*Y1( 1) 

DO 7 1=2,20 

YSTAR(I)=Y1( I)-(RH0HAT*Y1( 1-1)) 

7 CONTINUE 

C X TRANSFORM 

XSTAR( 1 , 1)=( l-( RH0HAT**2))**0. 5 
DO 9 J=2,K 

XSTAR(1,J)=((1-(RH0HAT**2))**0. 5)*X(1,J) 
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9 CONTINUE 

DO 11 L=2 ,T 

XSTAR(L,l)=l-RHOHAT 

11 CONTINUE 

DO 12 1=2, T 
DO 13 J=2, K 

XSTAR( I , J)=X( I , J)-RHOHAT*X( 1-1 , J) 
13 CONTINUE 

12 CONTINUE 

RETURN 

END 
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MSEB 



C THIS PROGRAM IS DESIGNED TO CALCULATE THE MEAN SQUARE ERROR OF 
C THE BETA VECTOR 

DIMENSION B1 ( 5000 ) , B2 ( 5000 ) , B3( 5000 ) , B4( 5000 ) , BS( 5000 ) ,B6(5000) , 
*B7( 5000), B8(5000),B9( 5000), B10( 5000), BX( 5000) ,BY(5000) 

OPEN ( UN IT=2 1 , F I LE= 1 DAT1 1 ) 

OPEN (UNIT=31,FILE='DAT2') 

OPEN (UNIT=41 , FI LE= ' DAT3 1 ) 

OPEN (UNIT=51 , FI LE= 1 DAT 4 1 ) 

OPEN (UNIT=61 , FI LE= ' DAT5 ' ) 

C 

COUNT=1000 

READ(21 ,900)(B1( I) ,B2( I) , 1=1,1000) 

CALL MSEBET (B1 ,B2, COUNT, XMSEDW) 

READ( 3 1 , 900 ) ( B3( I ) ,B4( I) , 1=1,1000) 

CALL MSEBET (B3,B4, COUNT, XMSETN) 

READ(41 ,900)(B5( I) ,B6( I) , 1=1,1000) 

CALL MSEBET (B5 ,B6 .COUNT, XMSEPW) 

READ( 51 ,900)(B7( I) ,B8( I ) , 1=1,1000) 

CALL MSEBET (B7,B8, COUNT, XMSEBM) 

READ( 61 , 900 ) ( B9( I ) ,B10( I) , 1=1,1000) 

CALL MSEBET (B9, BIO, COUNT, XMSEOLS) 

900 FORMAT (Fll. 8,2X,F11. 8) 

WRITE(6,*) 'MSEDW' 

WRITE(6 ,*)XMSEDW 
C 

WRITE(6,*) 1 MSETN 1 
WRITE(6,*)XMSETN 
C 

WRITE(6,*) 'MSEPW 1 
WRITE(6,*)XMSEPW 
C 

WRITE(6,*) 1 MSEBM 1 
WRITE(6,*)XMSEBM 



C 
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WRITE(6,*) 1 MSELS 1 

WRITE(6,*)XMSELS 

STOP 

END 

q **************** SUBROUTINE MSEBET **************************** 

c 

SUBROUTINE MSEBET(BX ,BY ,AN,XMSEB) 

DIMENSION BX( 5000) , BY ( 5000 ) , SUM( 5000 ) 

PLACE=0 
DO 901 1=1, AN 

SUM( I)=((BX( I)-1)*(BY(I)-1))**2 
P LAC E= PLACE + SUM(I) 

901 CONTINUE 

XMSEB= PLACE/AN 

RETURN 

END 
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